Equation of state of low— density neutron matter and the 1 S'o pairing gap. 
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We report results of the equation of state of neutron matter in the low-density regime, where 
the Fermi wave vector ranges from 0.4 fm -1 < fop < 1.0 fm _1 . Neutron matter in this regime is 
superfluid because of the strong and attractive interaction in the 1 So channel. The properties of 
this superfluid matter are calculated starting from a realistic Hamiltonian that contains modern 
two- and three-body interactions. The ground state energy and the 1 So superfluid energy gap are 
calculated using the Auxiliary Field Diffusion Monte Carlo method. We study the structure of the 
ground state by looking at pair distribution functions as well as the Cooper-pair wave function used 
in the calculations. 



I. INTRODUCTION 

Pure neutron matter is the natural first approxima- 
tion to the baryonic matter that composes the bulk of 
neutron stars. At very low densities below neutron drip, 
(i.e. where the Fermi wave vector is roughly, kp < 0.2 
fm -1 ) neutron star matter is conjectured to be nu- 
clei surrounded by a relativistic gas of electrons At 
higher densities the matter becomes liquid and very neu- 
tron rich. Here we study matter at Fermi wave vector 
0.4 fm _1 < kp < 1.0 fm -1 where it is reasonable to 
approximate it as pure neutron matter, and also extend 
some of our results into the lower density regime in order 
to compare with other calculations. 

At these densities the interaction is dominated by the 
1 5'o channel with a large and negative scattering length, 
a ~ —18.5 fm. The product of the effective range and the 
Fermi wave vector is of order unity, so, while the form of 
interaction cannot be neglected, it becomes less impor- 
tant. Analysis of the phase shifts of the neutron-neutron 
So interaction indicates that neutrons should pair and 
form a superfluid. Therefore the superfluid phase must 
be included when investigating the equation of state in 
this regime. 

Many methods have been used to approximately cal- 
culate the equation of state. One class of methods uses 
Skyrme or relativistic mean-field methods that use effec- 
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tive interactions that have been fit to the properties of 
nuclei. However, even those calculations that describe 
neutron-rich nuclei reasonably well, give rather different 
equations of state for pure neutron matter Q. We in- 
stead use a nonrelativistic Hamiltonian with two- and 
three-body interactions. All modern accurate two-body 
interactions fit the Nijmegen data[3| within experimen- 
tal errors and should give essentially the same equation of 
state at low density. At longer range, these interactions 
are dominated by the one-pion exchange and have strong 
spin-isospin dependence, which must be included for ac- 
curate predictions. Three- and higher-body interactions 
are less well known but in this density regime they are 
small. 

Our calculations extend the work we first reported in 
ref. The ground state of neutron matter is com- 

puted using the auxiliary field diffusion Monte Carlo 
Q (AFDMC) algorithm; it is an extension of the diffu- 
sion Monte Carlo methodp, and Green's function Monte 
Carlo method These Monte Carlo algorithms are very 
well suited to project a trial wave function onto the 
ground state in order to study the ground state prop- 
erties of a system. The Green's function Monte Carlo 
method has been used to study the properties of light 
nuclei with very high accuracy [g. The advantage of the 
auxiliary field diffusion Monte Carlo method over the 
Green's function Monte Carlo method is that it can be 
extended to larger nuclear systems; in fact it has been 
used to calculate properties of heavy nuclei [9|], neutron- 
rich is otop es fiol [ll[ and neutron [l^, G3 an d nuclear 
matter[14| by simulating systems with upwards of one 
hundred nucleons. 
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The equation of state of neutron matter in the low- 
density regime has been a subject of many previous 
calculations!, El M, EE E M, S M, 0, M m. 
While in this regime, different Hamiltonians and differ- 
ent methods give similar behavior for the energy as a 
function of the density, there are appreciable differences 
in other important properties. In particular, the value 
of the 1 Sq superfluid energy gap is at present not well 
clarified and strongly depends on the Hamiltonian and 
the solution methodjj]. In this paper we focus on both 
the energy and the energy gap by considering a fully re- 
alistic Hamiltonian, and solve for the ground state using 
the AFDMC technique. As a starting point for the cal- 
culation we considered two forms for the trial wave func- 
tion. The first is a filled Fermi sea having the properties 
of a normal Fermi liquid which we will call the normal 
phase. The second has neutrons paired in the 1 S'o chan- 
nel with a Bardeen-Cooper-Schrieffcr (BCS) superfluid 
structure [5^]. At a fixed density, we find that the su- 
perfluid phase of the system is only marginally favored 
compared to the normal phase. However, to calculate the 
superfluid energy gap, the BCS structure must be used. 



II. HAMILTONIAN 

We study the ground state of neutron matter beginning 
with the non relativistic nuclear Hamiltonian 
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where m is the mass of the neutron, and wy and Vijk 
are two- and three-body potentials. Such a form for the 
Hamiltonian (with the kinetic energy modified to take 
into account the mass difference of the neutron and pro- 
ton) has been shown to describe properties of light nuclei 
in a good agreement with experimental data (see ref. @ 
and references therein). All the degrees of freedom re- 
sponsible for the interaction between nucleons (such as 
the 7r, p, A, etc.) are integrated out and included in v.y 
and Vijk- 

At present, several realistic two-nucleon interactions fit 
scattering data with very high precision. We use the two- 
nucleon potentials belonging to the Argonne family [26| . 
Such interactions are written as 



$> p (r y )0W(z,j), 



(2) 



where 0"(i,j) are spin-isospin dependent operators. 
The number of operators M characterizes the interac- 
tion; the most accurate for the Argonne family is the 
Argonne AV18 with M=18 [26[. Here we consider a sim- 
pler form derived from AVI8, namely the AV8' [l?) with 
a smaller number of operators. For many systems, the 
difference between this simpler form and the full AV18 
potential can be computed perturbatively [H, [28[ , as has 



been done in all Green's function Monte Carlo calcula- 
tions to date. Most of the contribution of the two-nucleon 
interaction is due to one-pion exchange between nucleons, 
but the effect of other mesons exchanges as well as some 
phcnomcnological terms are also included. 

The eight 0^(i,j) operators in AV8' are given by the 
four central components I , t, • t j , er^ • (Tj , (tr^ • er j ) (t i ■ t j ) , 
the tensor Sy, the tensor-r component SV/T,; • Tj, where 
Sij = 3(<r,j ■rij)((Tj -Vij) —<Ji -a-j, the spin-orbit -Sij 
and the spin-orbit-r (Lij ■ Sij)(Ti ■ Tj), where and 
Sij are the relative angular momentum and the total spin 
of the pair ij. All the parameters describing the radial 
functions of each operator in AV18 are fit to nucleon- 
nucleon scattering data below 350 MeV in the Nijmegen 
database [|[. The AV8' interaction is obtained by start- 
ing from AVI8 and making an isoscalar projection. It 
is refit in order to keep the most important features of 
AVI8 in the scattering data and the properties of the 
deuteron \2l\ . 

The three-nucleon interaction is essential to overcome 
the underbinding of nuclei with more than two nucle- 
ons. While the two-nucleon interaction is fit to scattering 
data and correctly gives the deuteron binding energy, it 
is not sufficient to describe the ground state of light nu- 
clei with three or more nucleons. The Urbana-IX (UIX) 
potential corrects this, and was fit to obtain the correct 
triton energy using Green's function Monte Carlo and to 
correctly reproduce the expected saturation energy of nu- 
clear matter within the Fermi hypernetted-chain approx- 
imation [29l | . It contains a Fujita-Miyazawa term (30j 
that describes the exchange of two pions between three 
nucleons, with the creation of an intermediate excited A 
state. Again, a phenomenological part is added to sum 
all the other neglected terms. The generic form of UIX 
is: 

V ijk = V** + V R . (3) 
The Fujita-Miyazawaterm [30| is spin-isospin dependent: 

VjjTT = M-K ^ {-V. ,. A' ,;. | {t, • T^Tj ■ T k } + 
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-[Xij,Xj k ][Ti -Tj,Tj -T fe ] 



(4) 



where the Xij operators describe the one pion exchange 
(see ref. (3f| for details). The phenomenological part of 
UIX is 



Vifk = UoY.T^m^T^m^k) . 



(5) 



The factors A2-n and Uq are kept as fitting parameters. 
Other forms of three-nucleon interaction, called the Illi- 
nois forces [3l| . which includes three-nucleon Feynman 
diagrams with two-A intermediate states, arc available. 
Unfortunately they provide unrealistic overbinding of 
neutron systems when the density increases [l2l [32| a- n d 
they do not seem to describe, realistically, higher density 
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(i.e. p > po = 0.16 fm -3 ) nucleonic systems. However 
in the low-density regime considered in this paper, the 
contribution of the three-body interaction is very small 
compared to the total energy of the system, so that the 
small errors in the UIX interaction should have negligible 
contributions to the equation of state and energy gap. 



III. AFDMC METHOD AND THE PFAFFIAN 
WAVE FUNCTION 

Uniform neutron matter is simulated by solving the 
ground state of a fixed number N of neutrons in a pe- 
riodic box, whose volume is fixed by the density of the 
system. The ground state of the system is calculated by 
means of the AFDMC algorithm [jj. Diffusion Monte 
Carlo projects out the lowest-energy state from a trial 
wave function ipx by a propagation in imaginary time: 



</>(t) 
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where Et is a normalization factor. In the r — > oo limit 
the only component of tpT that survives is the lowest- 
energy one not orthogonal to ipT'- 



<p a = lim ip(r) . 



(7) 



The evolution in imaginary time is performed by solving 
the integral equation 



i>(R,r)= dR'G(R,R',T)<if) T (R'), 



(8) 



where G(R, R' ,r) is the Green's function of the Hamil- 
tonian that contains a diffusion term, coming from the 
kinetic operator in H , and a branching term from the po- 
tential. The exact form of G(R, R' , r) is unknown, but it 
can be accurately approximated in the limit of At — > 0. 
The above integral equation is then solved iteratively, 
with a small time step, for a sufficiently large number of 
steps. A detailed description of the algorithm as well as 
the importance sampling technique used to reduce the 
variance can be found in [33, [34j . 

The presence of spin operators in the Hamiltonian re- 
quires a summation of all possible good spin states in 
the wave function [35| . This summation grows exponen- 
tially with the number of neutrons; for example, for a 
system of 14 neutrons the computation of (ip(R)\ip(R)} 
is a sum of squares of 2 spin amplitudes. The explicit 
summation of spin states is performed in Green's function 
Monte Carlo, but not in AFDMC, calculations where the 
spin states are sampled using Monte Carlo techniques [j| . 
This sampling is performed by reducing the quadratic de- 
pendence of spin operators in the exponential to a linear 
form by means of the Hubbard-Stratonovich transforma- 
tion. The effect of an exponential of a linear combination 
of spin operators consists of a rotation of the spinor for 
each neutron during the propagation. In order to have 
an efficient algorithm, the trial function must be chosen 



so that it can be efficiently evaluated when each neutron 
is in a specific position and spinor state. 

Since both positions and spins can be sampled, the 
AFDMC method can be used to solve for the ground 
state of much larger systems - more than one-hundred 
neutrons - than Green's function Monte Carlo with full 
spin summations. 

More detailed explanations of the AFDMC method 
and how to include the full two- and three-nucleon in- 
teractions in the propagator can be found in refs. (l2l . 
Il3l [32l . I3II ] , where the fixed-phase approximation used to 
control the fermion sign problem is also discussed. 

The AFDMC method projects out the lowest energy 
state with the same symmetry as the trial wave function 
from which the projection is started. The general form 
of the trial wave function is 



§{R,S), 



(9) 



where R = (n, . . . ,rjv) represents the spatial coordi- 
nates and S = (si, . . . , sjy) the spin states of the neu- 
trons. The spin assignments Sj consist of giving the 
two spinor components for each neutron, namely the two 
complex numbers aj, bi where 



(10) 



and the {||), is the spin- up and spin-down basis. The 
function fj entering in the so called Jastrow part of the 
trial wave function has only the role of reducing the over- 
lap of neutrons and thereby reducing the energy variance. 
Since it does not change the phase of the wave function, it 
does not influence the computed energy value in projec- 
tions methods. The function fj is computed as described 
in ref. [3. 

The antisymmetric part $ of the trial wave function 
is usually given by the ground state of non-interacting 
Fermions (Fermi gas), which is written as a Slater de- 
terminant of single particle functions. For example ho- 
mogeneous systems are usually simulated by considering 
plane-waves as orbitals. In this case 

$„(i?, S) = A [<h.(ri>*i) ■ ■ ■ Mr N , sjv)] , (11) 
where A is the antisymmetrizer (see Eq. IA3|) . 

4> a {ri,Si) = e ika ^( Si \xs, ms ,a) , (12) 

and a is the set of quantum numbers of single-particle 
orbitals that are plane waves fitting the box. The correct 
symmetry of the ground state is given using the closed 
shells occurring when the total number of Fermions in a 
particular spin configuration is 1, 7, 19, 27, 33,... 

However, in superfluid neutron matter there is a strong 
coupling between neutrons, and a wave function heaving 
a BCS structure must be used. 

BCS pairing correlations can substantially change the 
nodal structure of a trial wave function 1371. 13811. This 
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change, which gives the off-diagonal long-range order of 
the superfluid phase, will greatly alter the fixed-phase (or 
constrained path) energy. In order to correctly describe 
the superfluid ground state with these quantum Monte 
Carlo methods, we need to use a trial wave function with 
explicit pairing. For central potentials and singlet pair- 
ing, the BCS trial function can be written as a determi- 
nant (37l. |39| . However, for problems with a tensor force, 
or for spin triplet pairing, a general pairing state must 
be used. 

A fully paired state of N neutrons can be written, as 
shown in appendix |A1 as 

■A [012034 ••• <Pn-1,n] , (13) 



Similarly, we can construct a general state with n paired 
and o unpaired orbitals for a total of N — 2n + o particles 
as 



A [</>12</>34 • ■ ■ 02n-l,2ri ■ ■ ■ 1pl(2n + 1) . . . l/j (N)] , (14) 



which is the pfaffian of the (N + o) x (N + o) skew- 
symmetric matrix [39j 
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where the lower o x o section is all zeroes. 
The pfaffian is the antisymmetric product 



have the form 



PL4 = A[a 12 a 34: a 56 ...a N - 1 



(16) 



The result is normalized such that every equivalent term 
occurs only once, and a,j = — a^. 

Just as the determinant of a dense matrix can be calcu- 
lated efficiently in order N 3 operations, similar elimina- 
tion methods can compute the pfaffian. The basic pfaf- 
fian calculational methods we use here, and have used, 
for all previous superfluid neutron matter studi es 3 , l40| , 
are described in some detail in section II of [4l|], and 
those results are summarized in appendix [5] along with 
some additional techniques needed for these nuclear cal- 
culations. 

The nuclear Hamiltonian has spin-dependent terms 
that can flip the spin. For the simpler case of a purely 
central potential, the Hamiltonian will not change the 
particles' spin. Therefore in this simpler case we can solve 
for the ground state in one sector where each particle has 
a specified spin, and we only need to antisymmetrize over 
the particles with the same spin. In that case, $bcs re- 
duces to a determinant. Since in our AFDMC method, 
the Hamiltonian can change the particles' spin, and the 
particles can then take on any spinor value, we need to 
be able to evaluate the trial wave function for arbitrary 
spinor values for each particle. Therefore the pfaffian 
which gives the full antisymmetric form must be used. 
As shown in appendix the pairing orbitals 4> we used 



Cq6 



(17) 



where the sum over a indicates the fc-space shells of the 
cube with k values 
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(n x x + n y y + n z z) 



(18) 



for integer n x , n y , and n z . The function \ is the spin- 
singlet wave function for two neutrons 

X ( Si ,s j ) = -j=(( Si s j \n)-(siSj\m) . (19) 

With the spin states given as spinors as in Eq. [10] this 
becomes 



b*a* 



V2 



(20) 



Note that if the pairing coefficients c a are zero for all 
\k a \ > kF, the pfaffian of Eq. [T5]is exactly the Slater de- 
terminant of spin up and down neutrons filling the Fermi 
sea, and the pfaffian form goes over to the normal liquid 
state. The parameters c a are chosen variationally by per- 
forming a correlated basis function calculation [40, [42| . 
However, various other wave functions were considered to 
ascertain the effect of a particular choice on the results. 
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kp [far 1 ] 


P [fm- 3 ] 


En/N 


Ebcs/N 


0.4 


0.00216 


1.289(2) 


1.239(2) 


0.6 


0.00730 


2.606(4) 


2.579(2) 


0.8 


0.01729 


4.277(7) 


4.305(3) 


1.0 


0.03377 


6.197(2) 


6.231(3) 



TABLE I: AFDMC energies per particle for 66 neutrons inter- 
acting with the AV8'+UIX interaction in a periodic box as a 
function of the Fermi wave vector and corresponding density 
p. The values E n correspond to the simulation of neutron 
matter using the Fermi gas ground state in the trial wave 
function, while Ebcs are the results obtained using &bcs- 
All the energies are expressed in MeV. 
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A. Equation of State 

We computed the energy of neutron matter by simulat- 
ing neutrons in a periodic box at densities corresponding 
to kp = 0.4, 0.6, 0.8 and 1.0 fm -1 using in the trial wave 
function both $ n and &bcs- We found that the abso- 
lute energy is slightly different depending on the choice 
of the trial function $. The results obtained using the 
two different trial wave functions are reported in table [J 
As can be seen, the BCS state is favored at &f = 0.4 and 
0.6 fin. - , while the normal state trial function gives the 
lowest energy at kp = 0.8 and 1.0 fm -1 . The maximum 
difference between the results for the two different trial 
wave functions is about 4 percent of the total energy at 
kp = 0.4 fm , probably because at such a low density 
the pairing between neutrons in the Sq channel is very 
important and <&bcs includes such correlations in the 
wave function in a more effective way. In the other cases 
the energies obtained with <&bcs arid &n are within 1 
percent. 

Since the coefficients entering in $ bcs were chosen by 
a correlated basis function calculation that adds a two 
body correlation factor to the usual BCS state [ifl |42|. 
in order to determine if this method is adequate for find- 
ing a good BCS form, we repeated some of the calcula- 
tions using different coefficients. In particular we tried 
using, as a pairing function, the solution from the uncor- 
rected BCS equation, as well as a pairing function with 
the same form as that of ref . [13, 0] • This calculation 
has carefully optimized coefficients, but the interaction 
is the 1 Sq channel of AV18 acting only between unlike 
spins. In both cases we find the energy is slightly higher 
than that found when using the correlated basis function 
coefficients. 

The equation of state of low-density neutron matter, 
computed using $bcs is displayed in Fig. [TJ compared 
to the diffusion Monte Carlo results of Gezerlis and Carl- 
son [22l | , to the variational cluster summation calculation 
of Friedman and Pandharipande [15( and to the results 
of Epelbaum, Krebs, Lee and Meifiner 2l|. The differ- 



FIG. 1: (color online) The equation of state of neutron mat- 
ter as a function of the Fermi wave vector kp. The energy 
has been divided by the energy of the noninteracting Fermi 

E-N. The AFDMC result is obtained us- 



3 ft k 



2m 



gas, E FG 

ing the full Hamiltonian AV8'+UIX (green circles), and com- 
pared with the results of Gezerlis and Carlson (red squares) 
who considered a simpler Hamiltonian [221 ] . The blue triangles 
correspond to the calculation of Friedman and Pandharipande 
using the Urbana vn two-nucleon interaction modified to in- 
clude some three-body effects [l5|]. The black diamonds show 
the results of Epelbaum, Krebs, Lee and Meifiner [2lj. In the 
inner part of the figure, the AFDMC energy in MeV is shown 
as a function of the density p in fm -3 is also displayed along 
with a curve to guide the eye. 



ences between the various calculations are due to differ- 
ent approximations and interactions used. The AFDMC 
method uses a realistic Hamiltonian containing a modern 
two-body and the corresponding three-body force. The 
variational cluster summation calculation was performed 
using the older Urbana U14 two-nucleon interaction[3| 
modified to include a density dependent term that mod- 
els the effect of a three-body force. As mentioned above, 
the calculation of Gezerlis and Carlson uses only the 1 So 
channel interaction of AV18 between unlike spins. This 
choice is motivated by the fact that this channel is dom- 
inant in neutron matter in this regime. However the ef- 
fect of other channels as well as using the 1 <So interaction 
partly in the triplet channel, since all unlike spin pairs 
interact, could play an important role in the many-body 
correlations of the system. Finally, Epelbaum and collab- 
orators computed the equation of state within the chiral 
effective field theory by simulating neutrons on the lattice 
up to the N 3 LO order pa}. 

Each of these calculations used different methods to 
solve for the ground state. Both the AFDMC and the 
diffusion Monte Carlo method used by Gezerlis and Carl- 
son are projection methods that, apart from the con- 
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straint used to control the Fermion sign problem, are ex- 
act. However, the constraint plays an important role in 
finding the correct ground state, and different trial func- 
tions give different constraints and therefore different re- 
sults. For these two Monte Carlo methods, trial functions 
have the same kind of BCS form. However Gezerlis and 
Carlson use a different approach for the choice of the co- 
efficients entering in the pairing orbitals of Eq. 1171 Their 
c a parameters are chosen by varying them to minimize 
the fixed-node energy [12, [13, Unfortunately this 

same technique is not currently applicable to AFDMC 
because the variance of the calculation is too high to be 
able to choose the coefficients in a reasonable amount of 
computational time. The c a used in our AFDMC calcu- 
lations arc chosen instead by using a correlated BCS wave 
function solved within the CBF/BCS theory [U SI as 
discussed above. The variational cluster summation cal- 
culation may suffer from important uncontrolled approxi- 
mations coming from the cluster expansion as we recently 
pointed out in our paper comparing the equation of state 
of neutron matter at higher densities [l3j |. In addition, 
the variational cluster summation calculation does not 
include any pairing correlations in the variational wave 
function. The equation of state can be computed using 
N3LO as described in [H, [2(| HH, and the results are 
available for a small number of neutrons (N=12). They 
predict an equation of state that is globally lower than 
the other results. This model, while very promising be- 
cause it attacks the problem from a more fundamental 
point of view, will need to be extended to larger systems. 

B. Superfluid gap 

In a full many-body calculation the superfluid gap can 
be evaluated by using the difference: 

A(N)=E(N)-±[E{N+l) + E(N-l)] , (21) 

where the number of neutrons N is taken to be odd. The 
AFDMC algorithm can be used to simulate very large 
systems with up to a hundred nucleons [12|, Gja, LLJ, [22] • 
Unfortunately, because the gap has to be evaluated as the 
difference between total energies of different systems, the 
statistical error related to A is proportional to the num- 
ber of neutrons, and we have not been able to develop 
an efficient method of correlated sampling. As a conse- 
quence, in principle, the number of neutrons is arbitrary 
but if N is too large, the statistical error affecting the 
gap becomes larger than the gap itself. The maximum 
number of neutrons used in this work is 68. 

Particular care was taken to check that the AFDMC 
had converged. The simulations were repeated with dif- 
ferent time steps. Neither the energy nor the gap is de- 
pendent on the time step used - the extrapolation to the 
zero limit is within our error bars. 

The gap is strongly dependent on the number of neu- 
trons for small N. For both kp = 0.4 and 0.6 fm -1 the 



A computed with N = 12 ... 18 is noticeably larger com- 
pared to that computed with N = 62 . . . 68. We find 
that at k F = 0.4 fm" 1 the gap is A(14) = 1.79(6) 
MeV and A(66) = 1.5(2) MeV, while at k F = 0.6 fin" 1 
A(14) = 2.59(6) MeV and A(66) = 2.1(2) MeV. This 
behavior is well described by the analysis of Gezerlis and 
Carlson who solved the BCS equation in the simulation 
cell, and then reproduced this effect by using diffusion 
Monte Carlo [12]. In their paper they calculate with up 
to 90 particles without observing a substantial change in 
the gap compared to that given by simulating the sys- 
tem with about 66 particles, giving us confidence that 
our gaps have converged. 



P [fin" 3 ] 




k F [fm '] 



FIG. 2: (color online) The 1 5'o pairing gap of neutron matter 
as a function of the Fermi wave vector kF computed with dif- 
ferent methods. In the figure we display works of Wambach et 
al. [H|], Chen et al. [IfJ, Schulze et al. wfli, Schwenk et al. 
Cao et al. [49]. Gezerlis and Carlson [22t| and Margueron et 
al. 50]. All the results are compared with a BCS calculation 
(dashed line). 



We report in Fig. [5] the superfluid gap computed with 
AFDMC using N = 62 . . . 68, compared with other calcu- 
lations. It is clear that the different methods used to com- 
pute the pairing gap give different results. The mean- 
field BCS result is essentially unchanged when other re- 
alistic two-nucleon interactions are used [12, [EH, and 
reaches a maximum of about 3 MeV. This is because all 
the two-body interaction are fit by reproducing the S- 
and P-wave components from experimental data. How- 
ever a realistic study of the pairing gap must include the 
corrections due to the polarization effects given by the 
medium. The various results can be essentially divided 
in two different groups according to the different way 
used to include this effect: the many-body calculation 
using effective interactions based on Brueckner theory 
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or Hartree-Fock calculations, or the microscopic calcula- 
tions (Monte Carlo methods or correlated basis function 
theory) where the whole Hamiltonian describing the sys- 
tem is solved. The many-body effective-interaction cal- 
culations of Wambach et al. [43 1, C hen et al. [4y], Schulze 
et al. [47| and Schwenk et al. [48| predict a large reduc- 
tion compared to the BCS gap, with a maximum gap 
of about 1 MeV. The microscopic calculations based on 
correlated basis function theory or using quantum Monte 
Carlo techniques show a reduction of the gap compared 
to the BCS result particularly at high densities, where 
the maximum is about 2.1 MeV using AFDMC and 2.4 
MeV with correlated basis functions. The other available 
quantum Monte Carlo result by Gezerlis and Carlson 
was performed for smaller densities because it neglects 
several contributions from other channels of the interac- 
tion [22l | . The recent results provided by other many- 
body techniques using Bruckner Hartree-Fock and new 
effective interactions by Cao et al. [49j and Margueron 
et al. [5(| predict a superfluid gap closer to the AFDMC 
result. Their maximum value of A is about 1.7 MeV. In 
addition, the different methods predict different densities 
where the gap reaches the maximum value. 



C. Pair distribution functions and pairing orbitals 

Besides computing energies, the structure of So pair- 
ing can be investigated by a qualitative study of pair dis- 
tribution functions. If the pair energy is low enough that 
only the 1 S'o or the 1 Sq and 3 Pi channels are important, 
the interaction can be written as v c (rij) + v a (rij)cri ■ <jj. 
Even though we keep the full interaction, it is interest- 
ing to look at the two-body distributions that have this 
form. The corresponding pair distribution functions are 
defined by 



and 



9c{r) 



1 



2irr 2 pN 



E 

i<j 



and 



9<r(r) 



1 



2irr 2 pN 



E 

i<3 



(V#> 



(22) 



(23) 



where p is the density. pg c (r)d 3 r is the probability of 
finding a neutron in an infinitesimal volume d 3 r at a dis- 
tance r from another neutron, while pga-(r)d 3 r is -3 times 
the probability of finding a neutron such that the two are 
in a singlet state plus the probability of finding a neutron 
such that the two are in a triplet state. In the limit of 
large r, g c {r) — > 1, while g a (r) — ► 0. 

Since over j is 1 in triplet and -3 in singlet channels, we 
can write singlet and triplet pair distribution functions, 
gs(r), where S = for the singlet and S = 1 for the 
triplet, 



1 



9o(r) = -,[9c{r) - ga{r)} 



(24) 



91 W = -^>9c{r) +g a {r)\ 



(25) 



Because AFDMC, like diffusion Monte Carlo, most eas- 
ily calculates mixed estimates 



(0) M = 



we extrapolate these from the variational values 



{Oh 



bfr\o\*fo 

(iPt\4>t} 



(26) 



(27) 



as(0)~2(0) M -(0)v. 

The pair distribution functions computed with 
AFDMC are shown in Fig. [3] Closed symbols refer 
to <? c (f) at various densities, while open symbols rep- 
resent g<j(r). The calculations were performed at differ- 
ent Fermi wave vector; black circles represents the g(r) 
at &f = 0.4 fm -1 , blue squares Uf = 0.6 fm -1 , red dia- 
monds kp = 0.8 fm -1 and green triangles kp = 1.0 fm -1 . 
As it can be seen, the strong interaction in the 1 Sq chan- 
nel is evident in both the g c {f) and g a (r) which exhibit 
a peak at the same distance. The peak value of g a (r) is 
about -3 times that of g c {r), and the peaks increase as 
the density is lowered. 
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FIG. 3: (color online) Pair distribution functions g c (r) and 
gait) for neutron matter as defined in the text. The curves 
with closed symbols are the g c (r), those with open ones indi- 
cate g,r(r) corresponding to different Fermi wave vector. See 
the text for details. 



The strong So correlation is more evident using the 
singlet and triplet channel distribution functions, which 
we show in Fig. 0] Closed symbols represent the singlet 
state of the pair, while open ones the triplet state at 



various Fermi wave vectors as indicated in the legend. 
The singlet channel becomes very strong and dominant 
when the density decreases. 
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FIG. 4: (color online) Pair distribution functions <?o(r) and 
<7i(r) as defined in the text. Closed symbols represent the pair 
distribution function projected into the singlet spin channel, 
while open ones the triplet spin channel. See the text for 
details. 



We can compare these pair distribution functions with 
those of a noninteracting Fermi gas, 



9E G (r) = \[l + l\r) 



and 



9 [ G (r) = l[l-l\r)] , 
where l(r) is the Slater function defined as 



l(r) 



i sin(fc^ r) — kp r cos(kpr) 



(28) 



(29) 



(30) 



We report in Fig. GUgoM (black circles) and 51(7") (red 
squares), and the corresponding g$ G {r) (green dashed 
lines) and gf G (r) (blue dashed lines) of the Fermi gas at 
Fermi wave vector kp = 1.0 fm^ 1 and kp = 0.4 fm -1 . 
The triplet pair distribution function does not differ very 
much from the noninteracting case; it does have a small 
deviation at large distances for kp = 1.0 fin -1 . This 
means that quantum correlations, in this channel, in this 
density regime, are not too important. They become rel- 
atively more important at higher densities. The singlet 
pair distribution function, instead, is completely differ- 
ent than that of the noninteracting Fermi gas. However, 
at kp = 1.0 fm -1 the peak of go{r) is not so very far 
from the maximum value of go (r) at the origin, while 





4 5 6 
r [fm] 



FIG. 5: (color online) Pair distribution functions go(r) and 
gi(r), as defined in the text, at /cf = 1.0 fm -1 (upper panel) 
and kF — 0.4 fm^ 1 (lower panel). See the text for details. 



at kp — 0.4 fm -1 the strong peak of the singlet is far 
from the noninteracting case. The singlet pair distribu- 
tion function is also compared with corresponding varia- 
tional correlated basis function calculations (black solid 
lines) u sing either Fermi hyper-netted chain (FHNC) ap- 
proach [52| for the normal phase or CBF/BCS [42| for 
the superfluid phase. It is evident that the strong peak 
of the singlet is due to presents of the strong correlations 
in the system. 

We plot in Fig. [5] the spatial part of the pairing 
function used in &bcs at kp = 0.6 fm -1 , along the 
three spatial directions 100, 110 and 111 obtained by 
using the correlated basis function coefficients. These 
are compared with the simulation cell Slater functions 
^ccii = J2k k<k F e%h r ■ The functions corresponding to 
each direction end at L/2, L/y/2 and ^/3/4L, where L is 
the side of the simulation cell. 
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FIG. 6: (color online) The spatial functions used in the pair- 
ing orbitals at kp = 0.6 fm _1 . The solid (blue) line is the 
function obtained using the correlated basis function (CBF) 
coefficients while the dotted (red) line is the simulation cell 
Slater function. 
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APPENDIX A: BCS WAVE FUNCTION 
PROJECTED TO FIXED N 

The original BCS wave function was not an eigenstate 
of particle number - i.e. it explicitly broke gauge sym- 
metry. For spin-singlet paired fermions with the pairs 
having total momentum zero, the BCS form can be writ- 
ten as 



|BCS) oc J] [u k 



10) 



(Al) 



where the a£ is the fermion creation operator for a par- 
ticle in the k wave vector and spin projection s state, 
with anticommutation relations 



We have reported a detailed computation of the equa- 
tion of state of neutron matter in the low-density regime 
where the system is superfluid and neutrons pair in the 
1 Sq channel. The superfluid gap was also computed. 
The presence of spin-dependent interactions means that 
the wave function must be written as a pfafnan of two- 
neutrons pairing orbitals, and the definition and the com- 
putation of the pfafnan was also discussed. 

The use of a realistic nuclear Hamiltonian without us- 
ing any effective interaction combined with the use of a 
very accurate projection technique makes these results 
benchmark for other methods. Because of the constraint 
used to control the Fermion sign problem, the results 
could in principle depend on the importance function. 
We carefully verified the effect of the wave function with- 
out observing a particular bias due to the fixed-phase 
constraint used in the calculations. 

We compared the computed equation of state with 
other results, and we observed important deviations that 
could be due both to the model Hamiltonian and to the 
methods used to solve for the ground state. We found 
that the So pairing gap is only somewhat lower than 
that predicted by the simple BCS theory for densities 
corresponding to kp < 0.5 fm , but the polarization 
effects due to the bulk are very important at higher den- 
sities where a large suppression of the maximum value of 
the gap with respect to the BCS prediction was found. In 
particular, the maximum value of the gap is a bit larger 
with respect to other recent calculations, and much larger 
than other calculations based on effective interactions. 



Whs, a i> s >} = 5k,k'5 s ,s' ■ (A2) 

The Uk and Vk here are functions only of the magnitude, 
k = \k\, and this spatial symmetry along with the fermion 
antisymmetry guarantees only singlet pairs. 

For our Monte Carlo calculations, it is simpler to use 
the projection of this state onto a fixed number of parti- 
cles, N, in a periodic simulation cell of side L. We write 
the antisymmetric position- and spin-projected states as 

A\r u s 1 ,r 2l S2,---r N ,s N ) 

= j-p (-l) P \P(ri,si,r2,S2,---r N ,s N )) 

permutations P 

= * V>£(ri)V>+>2)...^ s + (rjv)|0) , (A3) 

where P represents the permutation of the particle labels, 
and (— f) p is I (-1) for even (odd) permutations. The 
position and momentum creation operators are related 
by 



ks ^3/2 



j-L/2 


,L/2 


,L/2 


/ dx 


/ d y . 




l-L/2 J 


-L/2 J 


-L/2 



dze lk - r ifj+(r). 



The standard BCS state is usually normalized by 
choosing \uk\ 2 + \vk\ 2 = 1 . Since we will be project- 
ing out the part with N particles, even if we start with 
a normalized state, the projected part will no longer be 
normalized. There is then no advantage to taking a nor- 
malized state, and instead we divide by each of the Uk- 
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If one or more are zero, it simply means that we should 
drop the 1 term for that k value since it is always filled. 
We therefore take 



The particle-projected BCS wave function is then 



ibcs) = n 



1 I Vk + + 



|0> . 



(A5) 



bcs(R,S) = (H,5|BCS) = - 7 =(0\^ SN (r N )^ SN _ 1 (r N ^)...^ sl (r 1 )l[ 



Vk_ 



"fcT -hi 



|0) 



(A6) 



r 



This is readily evaluated using Wick's theorem [53] to 
change from the given order to the normal order. Con- 
tracting tp s (r) and a^ s gives 



(A7) 



From Eq. IA61 we see that either both a^f arL d a -ki m a 
pair or neither must be contracted with ips( r ) operators 
to give a nonzero result. One particular contraction oc- 
curs when ip Sl (t*i) and ip S2 (r 2 ) contract with such a pair 
in fei, il) S3 {r^) and V's4( r 4) contract with another pair in 
fc 2 , etc. This gives a term 

^e ik ^- r ^(s lS 2\ n) — e ik2<r3 ~ ri) (s 3 s 4 \ T|>... 

^l^ N Mr*-i-r N ) {sN _ lSNl U) t (A8) 
U k N / 2 

where we drop an unimportant overall normalization fac- 
tor. Choosing different k terms to contract with corre- 
sponds to summing over all values of the fei, fc 2 , etc. 
with the constraint that no two of the k n values should 
be equal (anticommutating two pairs of operators does 
not change the sign). Choosing other contractions com- 
pletely antisymmetrizes this form, and we can then in- 
clude all terms in the k sums since these cancel when 
antisymmetrized. The result is 

&BCS = A[(f)(ri,Si,r 2 ,S 2 ) ...(f>(r N ^i,S N ^i,r N ,S N )] . 

(A9) 

where, for spin-singlet zero-momentum pairs, 

4>{ru si,r 2 , s 2 ) = ]T ^ e * k < r ^ [( Sl s 2 \ U»] . 
k Uk 

(A10) 

Since the many-body antisymmetrizer will interchange 
the particles in </>, we usually explicitly antisymmctrize 
4>. We then get, up to an unimportant normalization, 

<Kn, si,r 2 , S2 ) - V ^e ik -^~^ [{s x s 2 \ T4) - <.s lS2 | |t) 
k Uk 

which explicitly demonstrates the singlet pairing. For a 
very large simulation cell, the spatial function would be 



spherically symmetric and therefore an 5* state. For the 
typical sizes of our simulation cells, the function has the 
symmetry of the cube as seen in Fig. [5] Other possible 
fully paired states have different <j>{r\, s\,r 2 , s 2 ), but still 
have the general form of Eq. IA9I 

Often we want to investigate systems which are not 
fully paired. Obviously, if we have an odd number of 
particles, at least one must be unpaired. We include 
unpaired particles in specific states by multiplying the 
|BCS) state by a product of creation operators (or linear 
combinations of creation operators) for those states. The 
only change to the particle number projection described 
above is that these creation operators must be contracted 
with one of the ipsir) or the result will be zero. For n 
pairs and o occupied single particle states, we have 

$BCS = A [()>U<j>U-<hn-l,2n1pl(2n + l)...1p (N)} 

(A12) 

which is Eq. [Mj 



APPENDIX B: PFAFFIAN CALCULATIONS 

Here we give some details on how to calculate the pfaf- 
fian. Proofs of the statements are given in ref. [4l| . 
The pfaffian of a skew-symmetric matrix has the follow- 
ing three properties: 

a. Multiplying a row and the corresponding column by 
a constant is equivalent to multiplying the pfaffian 
by a constant. 

b. Interchanging two different rows and the corre- 
sponding columns changes the sign of the Pfaffian. 

c. A multiple of a row and corresponding column 
added to another row and corresponding column 
does not change the value of the pfaffian. 

In addition, the matrix must have even rank for the pfaf- 
fian to be nonzero. Using these properties, it is straight- 
forward to use, for example, Gauss elimination to reduce 
the skew-symmetric matrix to a block diagonal form with 
2x2 blocks, whose pfaffian is just the product of the 
nonzero elements in the first superdiagonal. A Fortran 
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fragment showing the algorithm without pivoting for a 
complex matrix a of even rank n, is 

p=(l. 0,0.0) 
do i=l,n,2 
do j=i+2,n 

fac=-a(i, j)/a(i,i+l) 

a(i+l :n, j)=a(i+l :n, j)+f ac*a(i+l :n,i+l) 
a(j ,i+l:n)=a(j , i+1 :n)+f ac*a(i+l , i+1 :n) 
enddo 

p=p*a(i,i+l) 
enddo 

As in standard Gauss elimination, we search the current 
row for a large pivot element, and pivot using property 
b to bring this onto the superdiagonal so that we don't 
divide by small numbers a(i,i+l). 

At the same time, we calculate the inverse of the ma- 
trix. 

When one particle changes position or spin (or for 
calculation of one-body properties like the gradient, ki- 
netic energy or expectation of a spin operator) the skew- 
symmetric matrix A has one row and the corresponding 
column changed. Writing the matrix B to be equal to 
A except for the row k with new elements B^j and the 
corresponding column elements, Cayley showed [54| 



Pf[B] =Pi[A]J2 B kjA 



-l 
jk 



(Bl) 



For efficient algorithms with spin-dependent potentials 
we want to be able to change two particles. A straight- 
forward implementation would first change one row of 
the matrix as above and calculate the new pfaffian, and 



update the inverse. Then change the corresponding col- 
umn to obtain the skew-symmetric matrix and its inverse 
(its determinant is the square of the pfaffian obtained be- 
fore). This will require 0{N 2 ) operations. For each of 
the N second particles we will require 0(N) operations 
to calculate the new pfaffian if the first column is differ- 
ent for each pair. Unfortunately the result is 0(N A ) to 
calculate pairwise potentials. 

However, for our case, the operation needed on a col- 
umn or row is independent of the other column or row 
(except for the common element). We can therefore 
imagine doing a single update for particle 1 and using 
this for all the terms where the pair contains particle 1. 
The common element does not require an update and can 
be done separately. 

It is most efficient to write this as a set of matrix mul- 
tiplies. We define the new column j of the matrix to be 
Cij , corresponding to a spin or derivative operator on 
particle j. Defining 

k 

Gij = ^2 ^im^mk^kj = G ml P m j = ~Gji , (B2) 

mk m 

we find that the ratio of the new to old pfaffians with the 
two rows and columns denoted by i and j changed is 

= Ajl [A2f w + G^] + PuPjj - PuPa , (B3) 

where A new is the yl-matrix with new rows and columns. 
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